function diff = Diffusion(kappa, u)
% diff = Diffusion(kappa, u)
% 
% Diffusion calculates "div ( kappa grad u)". This is done in a symmetric
% fashion so that no direction is prefered. See [1] for details.
%
% input:    - u     : 2D image with scaler values
%           - kappa : Diffusion Tensor (scalar function))
%
% output:   - diff  : div ( kappa grad u)
%
% [1] M. J. Ehrhardt and S. Arridge, 
%     Vector-Valued Image Processing by Parallel Level Sets, 
%     IEEE Transactions on Image Processing (to appear), 2013
%
% 2013-07-25 --------------------------------------------------------------
% 
% Matthias J. Ehrhardt
% CMIC, University College London, UK
% 
% matthias.ehrhardt.11@ucl.ac.uk
% 
% -------------------------------------------------------------------------

    U = Extent(u, 'constant', 1);
    K = Extent(kappa, 'constant', 1);
    
    [M, N] = size(u);
    
    I = (1:M) + 1;
    J = (1:N) + 1;
    
    partial_x_forw_u = U(I+1, J  ) - U(I  , J  );
    partial_x_back_u = U(I  , J  ) - U(I-1, J  );
    partial_y_forw_u = U(I  , J+1) - U(I  , J  );
    partial_y_back_u = U(I  , J  ) - U(I  , J-1);
     
    mean_x_forw_kappa = 0.5 .* (K(I+1, J  ) + K(I  , J  ));
    mean_x_back_kappa = 0.5 .* (K(I  , J  ) + K(I-1, J  ));
    mean_y_forw_kappa = 0.5 .* (K(I  , J+1) + K(I  , J  ));
    mean_y_back_kappa = 0.5 .* (K(I  , J  ) + K(I  , J-1));
    
    diff = mean_x_forw_kappa .* partial_x_forw_u - mean_x_back_kappa .* partial_x_back_u ...
        + mean_y_forw_kappa .* partial_y_forw_u - mean_y_back_kappa .* partial_y_back_u;
    
    diff = diff(:);
end